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Abstract. The Closest Point Method for solving partial differential equations (PDEs) posed on surfaces was recently 
introduced by Ruuth and Merriman [J. Comput. Phys. 2008] and successfully applied to a variety of surface PDEs. In 
this paper we study the theoretical foundations of this method. The main idea is that surface differentials of a surface 
function can be replaced with Cartesian differentials of its closest point extension, i.e., its composition with a closest point 
function. We introduce a general class of these closest point functions (a subset of differentiable retractions), show that 
these are exactly the functions necessary to satisfy the above idea, and give a geometric characterization this class. Finally, 
we construct some closest point functions and demonstrate their effectiveness numerically on surface PDEs. 
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1. Introduction. The Closest Point Method is a set of mathematical principles and associated 
numerical techniques for solving partial differential equations (PDEs) posed on surfaces. It is an 
embedding technique and is based on an implicit representation of the surface S. The original 
method [ ] uses a representation induced by Euclidean distance; specifically, for any point x in 
an embedding space M" containing S, a point ecp(x) £ S is known which is closest to x (hence 
the term "Closest Point Method"). The function ecp is the Euclidean closest point function. 

This kind of representation of a surface S will be generalized: we introduce a general class of 
closest point functions (and denote a member by "cp") within the set of differentiable retractions. 
The closest point extension of a surface function u is given by M(cp(x)). The common feature of all 
closest point functions cp is that the extension M(cp(x)) locally propagates data u perpendicularly 
off the surface S into the surrounding space R". In this way, closest point extensions lead to 
simplified derivative calculations in the embedding space — ^because M(cp(x)) does not vary in 
the direction normal to the surface. More specifically, we have (for every closest point function): 

Gradient Principle: intrinsic gradients of surface functions agree on the surface with Cartesian 
gradients of the closest point extension. 

Divergence Principle: surface divergence operators of surface vector fields agree on the surface 
with the Cartesian divergence operator applied to the closest point extension of those vector 
fields. 

These are stated more precisely and proven as Principles 3.4 and 3.5. Combinations of these 
two principles may be made, to encompass higher order differential operators for example the 
Laplace-Beltrami and surface biharmonic operators. These principles can then be used to replace 
the intrinsic spatial operators of surface PDEs with Cartesian derivatives in the embedded space. 

Numerical methods based on these principles are compelling because they can re-use simple 
numerical techniques on Cartesian grids such as finite difference methods and standard interpo- 
lation schemes [ ]. Other advantages include the wide variety of geometry that can be repre- 
sented, including both open and closed surfaces with or without orientation in general codimen- 
sion (e.g., filaments in 3D [ ] or a Klein bottle in 4D [ ]). In this way, the Closest Point Method 
has been successfully applied to a variety of time-dependent problems including in-surface ad- 
vection, diffusion, reaction-diffusion, and Hamilton-Jacobi equations [18], [13], where standard 
time integration schemes are used. It has been shown to achieve high-order accuracy [ , ]. It 
has also been used for time-independent problems such as eigenvalue problems on surfaces [ ]. 

The remainder of this paper unfolds as follows. In Section 2, we review both notation and a 
calculus on embedded surfaces which does not make use of parametrizations. In Section 3, we 
define our class of closest point functions within the class of differentiable retractions. The key 
property is that, for points on the surface, the Jacobian matrix of the closest point function is the 
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projection matrix onto the tangent space. We show this class of functions gives the desired sim- 
plified derivative evaluations and that it is the largest class of retractions that does so. We thus 
prove all of these functions induce the closest point principles of [ ] (in fact, a more general di- 
vergence operator is established). Finally, we give a geometric characterization of these functions 
(namely that the pre-image of the closest point function intersects the surface orthogonally). Sec- 
tion 4 discusses general diffusion operators: these can be treated with the principles established 
in Section 3 by using two extensions but there are also many interesting cases (including the 
Laplace-Beltrami operator) where one can simply use a single extension. Section 5 describes one 
possible construction method for non-Euclidean closest point functions which makes use of a 
multiple level-set description of S. This construction method can be realized numerically which 
we demonstrate in an example. Finally, in Section 6 we use the non-Euclidean closest point rep- 
resentation to solve an advection problem and a diffusion problem on a curve embedded in W'. 

2. Calculus on Surfaces without Parametrizations. In this section, we review notation and 
definitions to form a calculus on surfaces embedded in R" (see e.g., [9, 16, 17, 2, 5, 6]). 

2.1. Smooth Surfaces. Throughout this paper we consider smooth surfaces S of dimension 
dim S — k embedded in R" {k < n) which possess a tubular neighborhood [ ]. We refer to this 
tubular neighborhood as B(S), a band around S. 

In the case that S has a boundary dS 7^ 0, we identify S with its interior S = S \ dS. Moreover, 
we assume that S is orientable, i.e., if codim S = n — k ^ there are smooth vector fields Nj : S ^ 
]R", i = 1, . . . ,n — k which span the normal space J\fyS — span{]Vi(i/), . . . , N„_;(-(i/)} (see e.g., [i. 
Proposition 6.5.8]). The vectors Ni{i/) are not required to be pairwise orthogonal, but we require 
them to be normed |N/(i/)| = 1. Finally, we define the R"^*^ -matrix N(i/) := {Ni{y)\ . . . |N„_j-(y)) 
which contains all the normal vectors as columns. 

The tubular neighborhood assumption is sufficient for the existence of retractions (see [ 0]). 
This is important because the closest point functions defined in Section 3 are retractions. Note 
that every smooth surface embedded in R" without boundary has a tubular neighborhood by [10, 
Theorem 5.1]. Moreover, the orientability side condition is not restrictive since it will be satisfied 
locally when referring to sufficiently small subsets of the surface. 

The matrix P{y) denotes the orthogonal projector that projects onto the tangent space TyS of 
S at y G S. P as a function of y G S is a tensor field on the surface and can be written in terms of 
the normal vectors as 

PiS^R"""", P{y) = I-N{y)-N{y)^ 

where N(y)^ denotes the pseudo-inverse of N{y). If codim S = 1 then P is given by P{y) = 
I — N{y) ■ N(y)^ since the matrix N{y) has only a single normed column vector. 

2.2. Smooth Surface Functions and Smooth Extensions. Given a surface S, we consider 
smooth surface functions: a scalar function w : S — )• R, vector-valued function f : S ^ R'", and 
vector-field g : S — > R". We call g a vector-field because it maps to R", the embedding space of 
S, and hence we can define a divergence operation for it. The calculus without parametrizations 
for such surface functions is based on smooth extensions, defined in the following. 

Definition 2.1. (Extensions) We call ue ■ 'Ran extension of the surface function m : S ^■ 

R — and likewise for vector-valued surface functions — if the following properties hold: 

• Oe c R" is an open subset of the embedding space. 

• contains a surface patch S H Ci£ 7^ 0. 

• "cIsnriE = "Isnng/ the extension and the original function coincide on the surface patch. 

• M e C'(S n Of) ^ Me G C'(n£), the extension is as smooth as the original function. 
As extensions are not unique, ue denotes an arbitrary representative of this equivalence class. 

Different extensions might have different domains of definition. Here is just a generic 
name for an extension domain that suits the chosen extension and that contains the surface point 
y E S which is under consideration. 

An important point is that such extensions exist and may be obtained by using the following 
feature of a smooth surface: for every y E S there are open subsets n£,W of R", where y G Of, 
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and a diffeomorphism ip : W {xp,xp are as smooth as S) that locally flattens S 

(/^(ngnS) = WnlR^, fc = dimS, := {x e M" : X = (xi,...,x,„0,...,0)} , 

see e.g., [11, Chapter 3.5]. Now, : W n ]R§ ^ CIe n S parametrizes the patch n S of S, so 
u o ip~'^ : W n Rg ^> R is a smooth function. Let P E W^" be the projector onto Rq. For x e D,£, 
the function P ■ i/'(x) maps smoothly onto W n Rq, hence 

Me : Oe ^> R , me(x) := m o ip~^ (P ■ \p{x)) 

is smooth and extends u as desired. 

2.3. Calculus without Parametrizations. Now, we define the basic first order differential 
operators on S without the use of parametrizations (compare to e.g., [9, 16, 17, 2]): 

Definition 2.2. Let S be a smooth surface embedded in R", y e S a point on the surface, and 
P{y) G R"^" the projector (onto TyS) at this point y. Then the surface gradient Vg of a scalar function 
u, the surface Jacobian Dg of a vector-valued function f, and the surface divergence divg of a vector-field 
g are given by: 

Vsu{yY -VuT.iyY -Piy) , 

Dsfiy) := DfE{y) ■ P{y) , 
divs^(y) := trace(Dsg(y)) = trace(DgE(y) ■ P{y)) ■ 

Here, V fs the gradient and D the Jacobian in the embedding space R" applied to the extensions of the 
surface functions. The extensions are arbitrary representatives of the equivalence classes of Definition 2.1. 

Remark. There are other works (see e.g., [ ]) that use a variational definition of the surface 
divergence operator, which we denote by divg g, as 

j V dw*sgd'H"{y) = - J {Vsv,g) dn"iy) , e c^{n) , n c S . 
n n 

But this definition applies only to tangential vector-fields g, because this tangency is required by 
the surface Gauss-Green Theorem. The connection to our definition is as follows 

divsg = divs(Pg), 

because divg takes into accoimt only the tangential part Pg of the vector-field g. The two are equal 
if the vector-field g is indeed tangential to the surface. The trace-based definition of surface diver- 
gence in Definition 3.1 also applies to non- tangential fields [ ]. For example when codimS = 1 
and g is a non-tangential vector-field, it can be shown that 

divsg = divs(Pg) +divs((N,g) ■ N) = divs(Pg) + {N,g) ks , 

where N is the normal vector, and Kg = divg N is the mean curvature of S. That is, the surface 
divergence of g is the surface divergence of the tangential component plus an extra term which 
depends on the curvature(s) of S. 

3. Calculus on Surfaces with Closest Point Functions. Our aim is to compute surface in- 
trinsic derivatives by means of closest point functions. Closest point functions are retractions with 
the essential property that their Jacobian matrix evaluated on the surface is the projector P. 

Definition 3.1. (Closest Point Functions) We call a map B(S) S a closest point function cp, if 

1. cp is a differentiable retraction, i.e., cp :B{S) ^ S features the properties 

a) cp o cp = cp or equiv. cp | g = idg 

b) cp fs differentiable 

2. Dcp{y)=P{y)forallyeS. 
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If cp belongs to this class of closest point functions, we recover surface differential opera- 
tors (those in Definition 2.2) by standard Cartesian differential operators applied to the closest 
point extensions m o cp. This fundamental point is the basis for the Closest Point Method and is 
established by the following theorem: 

Theorem 3.2. (Closest Point Theorem) If cp : B(S) S is a closest point function according to 
Definition 3.1, then the following rule 

D[/ocp](y) = Ds/(y), y E S, (3.1) 

holds for the surface Jacobian matrix Dgf of a smooth vector-valued surface function f : S ^ R"'. 

We also show that any differentiable retraction that satisfies the relation in (3.1) must be a 
closest point function satisfying Definition 3.1. 

Theorem 3.3. Let r : B(S) S be a differentiable retraction (i.e., only part 1 of Definition 3.1 is 
required). If for every smooth surf ace function f : S ^ IR™ the following holds 

D\for]{y) = Dsfiy), y e S, (3.2) 

then r satisfies also part 2 of Definition 3.1 

Dr{y)=P{y) for all y E S, 

and is thus a closest point function. 

Proof. Ad Theorem 3.2: the closest point extension can be written in terms of an arbitrary 
extension /f: 

/ocp(x) =/£ocp(x) , X e B(S) . 
Now, we expand the differential using the chain rule^ 

D[focp]{x) = D/eocp(x) ■ Dcp(x) , X e B(S) . 

Finally, we set xtox = yESa point on the surface and read off the surface differential as of 
Definition 2.2. Because by Definition 3.1 we have cp(i/) = y and D cp(i/) = P{y), hence 

D[/o cp] (y) = D/E(y)-P(y) = Ds/(y) , yeS. 

Ad Theorem 3.3: we consider the identity map on S, idg : S ^ S, ids(y) = y. The sim- 
plest extension of idg is given by the identity map on the embedding space, id : R" — > R". By 
Definition 2.2, the surface Jacobian of the surface identity is 

Dsids(y) = Did(y)-P(y) =P(y), yeS. 

By assumption, (3.2) is true for every smooth surface function. And so with / = ids ^r^d the 
previous result we have 

Dr(y) = D[idsor](y) = Dsids(y) = P(y) , yeS. 

□ 

Direct consequences of the Closest Point Theorem 3.2 are the gradient and the divergence 
principles below. 

Principle 3.4. (Gradient Principle) If cp : B(S) — > S is a closest point function according to 
Definition 3.1, then 

V[Mocp](y) = VsM(y) , yeS, (3.3) 

holds for the surface gradient V su of a smooth scalar surface function m : S ^ R. 

Principle 3.5. (Divergence Principle) If cp : B(S) S is a closest point function according to 
Definition 3.1, then 

div[g o cp] (y) = div sg{y) , y e S, (3.4) 

holds for the surface divergence divgg of a smooth surface vector-field g : S ^ R". Notably the vector 
field g need not be tangential to S. 

^In this paper, we follow the convention that differential operators occur before composition in the order of opera- 
tions. For example, D/e o cp(x) means differentiate /g and then compose with cp(x). 
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3.1. Geometric Characterization of Closest Point Functions. A further characteristic feature 
of closest point functions is that the pre-image cp~^ (y) of a closest point function must intersect 
the surface S orthogonally. This fact will be established in Theorem 3.7, for which we will need 
the smoothness of the pre-image. 

Lemma 3.6. Let r : B(S) Sbea dijferentiable retraction and let y & S. The pre-image r~^{y) is 
(locally around y) a smooth manifold. 

Proof. By differentiation of the equation idg or{x) = r{x) and setting x = y G S thereafter, we 

get 

Por{x)-Dr{x) = Dr{x) =^ P{y) ■ Dr{y) = Dr{y) . (3.5) 

We see that the range of the Jacobian matrix Dr{y) is contained in the tangent space TyS. On the 
other hand, we can differentiate r o idg (y) — idg (y ) w.r.t. to y e S and right-multiply by a matrix 
T{y) G ]j^«x/c Yvhich consists of an orthonormal basis of 7^ S to get 

Dr{y)-Piy)=P{y) Dr{y) -Tiy) = T{y) . (3.6) 

Equation (3.5) shows that the linear operator Dr[y) : R" ^ R" is rank deficient while Dr{y) : 
TyS TyS is full rank by (3.6) and the rank is A: = dim S. Now, we can apply the Implicit 
Function Theorem to 

fy{x) = T{yf-{r{x)-y)=0 

fy'.W^ R*^, where r~^{y) is the set of solutions and Dfy{y) = T{yY ' ^^iy) is full rank, and 
hence r^^ (y) is smooth. □ 
Theorem 3.7. Let r : B{S) ^ S be a differentiable retraction and let y & S. The retraction r has the 
property Dr{y) = P{y) (and hence r is a closest point function according to Definition 3.1), if and only if 
for every y ^ S the smooth manifold r~^{y) intersects S orthogonally, i.e., 

Tyr-\y)=HyS. 



Proof. Let y & S and assume that the pre-image (y) intersects S orthogonally. Let ^ : 
(—£,£) — 7> r~^{y), e > 0, be a smooth regular curve in r~^{y) with ^(0) = y. As r~^{y) is smooth 
by Lemma 3.6 such curves ^ exist. Because ^(f) G r~^{y), r maps every point ^(f) to y: 

ro^[t)=y, for all f e (— £, e). 

Differentiation of the latter equation w.r.t. t and the substitution of f = yield 

Dr(y)-^'(0) =0. (3.7) 

By assumption ^'(0) is an arbitrary vector of the normal space MyS. Let (here) N(y) e '^nx(n-k) 
be a matrix where the columns form an orthonormal basis of MyS, then (3.7) implies 

Dr(y) ■ N(y) = . 

By the latter result combined with (3.6) the product of Dr{y) with the n x n orthogonal matrix 
(T(y) |N(y)) (the matrix T(y) E R"^*^ consists of an orthonormal basis of 7yS) is 

Dr(y)-(T(y)|N(y)) = (T(y)|0) ^ Dr{y) = T{y) ■ T{yf = Piy) . 

This proofs the first direction. 

Let now Dr(y) = P{y) for all y E S, i.e., the retraction r is already a closest point function. 
Let again f : (—£,£) — >• J'~^(y) be a smooth regular curve in r^^{y) with ^(0) = y. We consider 
again (3.7) but this time we replace Dr{y) with P{y): 

P(y)-^'(0)=0. 

The latter tells us that every vector of Tyr^^ (y) is perpendicular to TyS, in other words r~^ (y) 
intersects S orthogonally. □ 
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3.2. The Euclidean Closest Point Function. The class of closest point functions from Defini- 
tion 3.1 is not empty. We show that the Euclidean closest point function — which is well-defined 
by our assumptions on surfaces in Section 2.1 — belongs to this class. However, this is not just a 
corollary of the last theorem, we have to show that ecp is differentiable on the surface S. 

Theorem 3.8. The Euclidean closest point function ecp is a closest point function satisfying Defi- 
nition 3.1. 

Proof. The proof is instructive in the special case of codim S = 1 and we begin with that case: 
The Euclidean closest point function ecp : B(S) Sisa continuous retraction and can be written 
in terms of the Euclidean distance map d as 

ecp(x) = X — d{x) ■ Vrf(x) . (3.8) 

Because codim S — 1, we can replace d with a signed Euclidean distance sd which on B(S) is a 
classical solution of the Eikonal equation, where sd is as smooth as S and where V sd{y) = N{y), 
for y £ S, is exactly the surface normal. So ecp is differentiable with Jacobian matrix 

Decp(x) = I - Vsd(x) ■ V sd(x)^ - sd(x) ■ D^sd{x) , 

and, as sd vanishes on S, we get 

Decp{y) = I-N{y)-N{yf = Piy), yeS, 

which proves the statement in this case. 

If S is of higher codimension, we cannot use this argument since d is not differentiable on the 
surface S. But it is differentiable off the surface and so is ecp by (3.8). We prove now that D ecp 
extends continuously onto S and equals the projector P there. 

Let 7 : U C R*^ — )• Og, 6i 7(61) be a regular parametrization of a surface patch Qg C S, 
(k = dim S). Then, by adding the normal space 

O(0) =7(0i) + No7(0i).02, 

O : n c R" ^ B(S) , = {01,62) ^ ^{0) , 

we get a coordinate system on a corresponding subset of B(S). N{y) = {Ni{y)\ . . . \N„_i^{y)) is a 
matrix formed of the normal space basis and 62 E R"^*^. Now, we have 

ecp 00(0) = 7(0^) . 

As long as 02 7^ (i.e., off the surface) we can apply the chain rule: 

DecpoO- {Dg^y + D0^{N o j) ■ 02 | N07) = {Dg^y \ 0) . 

The the second factor of the left hand side extends onto S (i.e., 62 = is admissible) and is 
invertible, the right hand side is defined for 02 = anjrway. After inverting, we send 62 to zero 
and so we have (0(01,0) = 7(0i)) 

Deep 07 = {De^j \ 0) ■ {Dg^'y \Noj)-\ 

As Dgj 7 and N o 7 are orthogonal, i.e., Dg^ 7^-^07 = and N o 7^ • Dg^ 7 = 0, we can write the 
inverse in terms of the pseudo inverses of the sub-matrices 

and so 

Deep 07 = (Dej7 I 0) ■ {Dg^y \ Noj)-^ = {Dg^j \ 0) ■ (^^^ = Dg^j ■ Dg^ = P07 . 
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In the latter equality we have used that the outer product of a full column-rank matrix A and its 
pseudo inverse gives the projector onto the image of A which in our case is the tangent space. 
Finally we get the assertion: 

Decpo7 = Po7 D ecp(i/) = P(i/) , J/ e S . 

□ 

Remark. As a consequence of Theorem 3.8 we see that the function 

X — ecp(x) = d{x) ■ Vd(x) = V 

is differentiable. The article [8] introduces v{x) := x — ecp(x) as the vector distance function 
and refers to the article [ ] for some of its features. In particular the authors of [ ] prove that 
the squared distance function is differentiable and show that the Hessian matrix D'^rj[y) of the 
function y] := when evaluated at a surface point y G S is exactly the projector onto the 
normal space MyS. Here, we have an alternative proof that rj (and thereby the squared distance 
function) is differentiable with derivative I — P{y) = Dv{y) = D^r] (y) on S which is the projector 
onto MyS. 

4. Surface Intrinsic Diffusion Operators. In this section we discuss the treatment of diffu- 
sion operators, that is, second order differential operators of the form 

divs(A(y)Vsw), 

in terms of the closest point calculus. Of course, by combining Principles 3.4 and 3.5, we can set 
up the diffusion operator as follows 

w := A o cp V[m o cp] 
^ j^(y) = A(y)VsM(y) , yeS 
=> div[a;ocp](y) =diwsw{y) = divs(A(y)VsM(y)) . 

In this set-up, we have a second extension z<; o cp in the last step. liv = A(y)VsM is a tangential 
vector field, we call the operator a surface intrinsic diffusion operator (these are relevant regarding 
the physical modeling of surface processes). The subject of this section is that we can drop the 
second extension in many cases (depending on A and cp) given the tangency of v. The key to 
this result is the following lemma on the divergence of vector fields which are tangential on the 
surface while tangency is allowed to be mildly perturbed off the surface. 

Lemma 4.1. (Divergence of tangential fields) Let She a smooth surface, and let j : U C 'K'^ ^ 0.$ 
(k = dimS) be a regular parametrization of a generic surface patch Cig c S. Let CI be a corresponding 
open subset of the embedding space ]R" such that Q n S = Qg. Let ^ : U xV ^ Cl(V c IR"^*^) be a 
coordinate system on O that satisfies 

cD(0i,02) =7(0i) + No 7(0i) .02 + 0(1021^), as 02^0. 

Here the columns of the n x {n — k) matrix N{y) := (Ni(y)| . . . |N„_jt(!/))/ V & S, give some basis of 
the normal space MyS. Let v : B{S) ^ W be a vector field and let 

v{ei,02) = z;ocD(0i,02) = De^7{ei)-v{ei,92)+Noj{9^) • t]{9i,e2) 

be the representation of this vector field on Ciin 9 = {9i,92)-variables with a decomposition into the two 
components (v E IR*^, rj E which are tangential and normal to the surface S. 

If the coefficient r] of the normal part satisfies 




r]{9i,0) = and trace Dg^rj (61,0) = , 
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then (as the surface patch is generic) 

drvv{y) — divg v{y) , y & S . 

The condition r}{d\,Q) = means that the restriction v\s is a tangential surface vector field. 

Proof. We show divs vo'y — div vo'y. For a tangential field the surface divergence in terms 
of a parametiization 7 is given by 

divsi;o7 = J-dweM^i,Q)^) = dwe,{v{^i,Q)) + £ |f i?;(ei,0) . (4.1) 

where g = det G and G = Dg^ 7^05^7 is the metric tensor induced by 7. 
Next, we turn to div o 7. In a first step, we obtain 

divi;o7 = tiaceDi;o7 = tiace Di? o *(0i,O) = trace {pev{ei,0) ■ D0*(0i,O)"^) . 

The derivative of v is 

k n—k k n—k 

D0v{9i, 62) = £ dijVevJ +l^Nio jVetJi +l^vi diDej +l^tji De{Ni o 7) . 
/=i /=i ;=i /=i 

k n—k k 

De5(0i,O) = £ dijVevJ +^Nio jVevI + L ^' ^'^eT 
1=1 1=1 ;=i 

since j; {61, 0) = 0. By oiir assumptions on <E> we have 

De^{0i,O) = {De,j\Noj), De^{9,,or' = {De,7 \ N o j)-^ = 
Now, we write the divergence as 

k , , n—k 



divi;o7 = £^trace (dijWevJ ■ Dg^{0l,O)~'^^ + ^ trace (NiojWerjf ■ De^{0i,Oy^ 
/=l / , , V 

k . 

+ £trace(a/D07 ■ Dq<^{Qx,^)-^ 



k 

/=1 

Using the rules of the tiace operator we get 

k 



divi;o7 = trace [l)^ Q 0^0(01, O)-HD0j7 | No7)j +^^^trace (a;De7 • D0cD(0i,O)-i) i;. 
Next, we use D07 = (09^7 | 0) and the second assumption tiace DQ^r]{9\, 0) = to get 
divi7 o 7 = trace {dq ) + ^"^^^^ ' ^e-yl'^ ^\ 



;=i 

k 



: trace {Dq^v) + ^ trace (3/00^7 • 00^7+) i); (4.2) 

k . . 

: div0j(i;(0i,O)) + ^ trace (3/00^7 • Df,j7+j 0/ 



The final step is to show that tiace (3;Dej7 • 09^7''') = Partial differentiation of the definition 
of G = Dg^j^De^j and using the definition of 0$^^^ — G'^De^j^ yields 

trace (diDe.j ■ De.j^) = \ trace(G-ia,G) (4.3) 
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While partial differentiation of the equation G G = I results in the (matiix) ordinary differential 
equation (ODE) 

di (g-1) =-(G-i3,G) G-\ 

which unplies 3/ det G^^ = — trace(G^^3;G) det G~^ and hence g = det G satisfies 

-\dig = di^ =3;detG-i = -trace(G~ia;G)detG"^ = - trace(G~^3;G) ^ . 

This result combined with (4.3) completes the proof: we replace the last term in the last equality 
of (4.2) with 

trace (diDe.j ■ De.i') = ^ trace(G-ia,G) = ^ , 

and compare with (4.1) to see that divg c o 7 = div c o 7. □ 

4.1. Diffusion Operators with Less Extensions. Armed with Lemma 4.1 we obtain a closest 
point calculus involving less extensions, for all surface intrinsic diffusion operators. There are 
two situations depending on how much "help" we get from the diffusion tensor A. First, we 
look at the situation where maps all vectors ^ G IR" to tangent vectors. 

Theorem 4.2. Let A : S ^ ^nxn ^ tangential diffusion tensor field, that maps all vectors to 
tangent vectors, i.e.. My & S we have e M" ^ ^(y)? G TyS. Then the corresponding surface 
diffusion operator can be written in terms of a closest point extension as 

divs (AVgw) (y) = div (A o cp V[m o cp]) (y) , y & S . 

A re-extension 0/ V[m o cp] fs not necessary and the result is true for all closest point functions. 

Proof. As A maps all vectors to tangent vectors the field v := A o cp V [m o cp] is tangential to S 
for all X. Now, we want to apply Lemma 4.1, that means we have to define the coordinate system 
O and the field v{Qi,92). We reuse the parametrization 7 of a generic surface patch from Lemma 
4.1. For a fixed value of 61 and the corresponding surface point 7(61) we simply parametrize 
(with 62) the pre-image cp^^(7(0i)) in order to get ^{9i, ^2). Since by Theorem 3.7 the pre-image 
cp~^(7(0l)) intersects S orthogonally, we can organize its parametrization in such a way that O 
satisfies 

0(01,02) = 7(^1) 7(^1) ■ ^2 + C'(|02|^) , as 02 ^ . 

As our particular vector field v is tangential to S for all x in a band B(S) around S, the corre- 
sponding V from Lemma 4.1 takes the form 

v{er,e2) = vo^{ei,e2) = De,^{ei)-v{er,e2) 

with rj equal to zero. Hence, by Lemma 4.1 we conclude the first equality of 

div(A o cp V[m o cp])(y) = divs(A o cp V[m o cp])(y) = divs(AVsM)(y) , y e S , 
while the second equality is because we have cp I5 = ids and V[w o cp] |g = Vgw. □ 

Next, we look at the situation where A{y) maps only tangent vectors ^ e 7^S to tangent 
vectors. This case is particularly interesting since — by considering diffusion tensors of the form 

A{y)=a{y)-1, (4.4) 

where 1 G M" ^ " is the identity matrix on the embedding space — it covers the special case where 
we are given a scalar diffusion coefficient a{y) G R. It is clear that the diffusion tensor in (4.4) can 
only map tangent vectors to tangent vectors (this is weaker than the requirement in Theorem 4.2). 

Theorem 4.3. Let A : S ^ ^nxn ^ tangential diffusion tensor field, that maps only tangent 
vectors to tangent vectors, i.e., My & S we have e TyS =^ ^(j/)? G TyS. If cp is a closest point 
function such that the transpose of its Jacobian matrix maps tangent vectors to tangent vectors, i.e., 

Vx e B(S) wehave: e T^pj^-jS =^ D cp(x)"''^ e T^p^^jS . 
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then the corresponding surface diffusion operator can he written in terms of a closest point extension as 

divs (AVsm) (y) = div(A o cp V[m o cp])(i/) , y e S . 

A re-extension 0/ V [m o cp] fs not necessary. 

Proof. As A maps only tangent vectors to tangent vectors the field v := A o cp V[m o cp] will 
be tangential to S only if V[m o cp] is tangential to S. Since cp as a closest point function satisfies 
cp = cp o cp and D cp o cp = P o cp we get 

D cp = D [cp o cp] = D cp o cp D cp = P o cp D cp . (4.5) 

Now, we expand V[m o cp] by referring to an arbitrary extension Ue of u as of Definition 2.1: 

V[m o cp] = Dcp"'^ Vme o cp = Dcp^ P o cp Vm£ o cp = Dcp"'^ Vgw o cp . 

Since ^ : = Vg m o cp (x) belongs to T^p^x) S, V [m o cp] (x) also belongs to 7^p(x) S by our requirement 
onDcp(x)^. Now, we are sure that c := AocpV[Mocp] is tangential to S for all x, and so the 
same argumentation as in the proof of Theorem 4.2 applies. □ 

4.1.1. Remarks. 

• If cp is a closest point function with a symmetric Jacobian matrix D cp = D cp^, then D cp^ 
has the property required in Theorem 4.3: (4.5) combined with the symmetry yields 

P o cp D cp = D cp = D cp^ 

which shows that D cp^ in this case maps all vectors ^ to tangent vectors. An example of such a 
closest point function with a symmetric Jacobian matrix is the Euclidean closest point function 
ecp (see Theorem 3.8): 



Decp(x) 



I - Vd{x) ■ Vd{xY - d{x)D^d{x) x e B(S) \ S 
P(x) xeS 



• In both theorems the crucial bit is that A{y) for fixed y maps certain vectors ^ to the tangent 
space, i.e., A{y)^ e TyS. The theorems are still true if we let A also depend on the function u, 
for example this dependence could be of the form A(y, M(y), VsM(y)). 

4.2. The Laplace-Beltrami operator. As an example we discuss the surface Laplacian (the 
Laplace-Beltrami operator) in terms of the closest point calculus: 

1. We may write the Laplace-Beltrami operator as 

Agw = divs(Vsw) . 

Given a closest point function that satisfies the requirement of Theorem 4.3, e.g., the Euclidean 
closest point function ecp, we have 

A[Mocp](y) = AsM(y) , yeS. 

2. Alternatively, we may rewrite the Laplace-Beltrami operator by using a diffusion tensor as 

Asm = divs(PVsM) . 

The diffusion tensor here is the projector onto the tangent space. So we have replaced the iden- 
tity matrix with the surface intrinsic identity matrix. Now Theorem 4.2 allows us to compute 
the Laplace-Beltrami operator like 

div(P o cp V[m o cp])(y) = AsM(y) , yeS, 

without any further requirements on D cp. 
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3. If P is not known a priori and cp does not satisfy the requirement of Theorem 4.3, we can 
always work with re-extensions (directly applying Principles 3.4 and 3.5) 

div(V[M o cp] o cp)(y) = A5m(i/) , y e S. 

Note by expanding the expression 

V[m o cp] o cp = (Dcp"'' Vm o cp) o cp = (Dcp o cp)"'^VM o cp 
= (P o cp)^VM o cp = P o cp Vm o cp, 

we can see that re-extension in fact means an implicit version of approach 2. 

4.3. Proof of a Principle of Ruuth & Merriman. For the original Closest Point Method (with 
ecp) Ruuth & Merriman in [ ] also reasoned that a re-extension is not necessary for some diffu- 
sion operators based on the fact that the Euclidean closest point extension satisfies the PDE 

(Vd(x),V[Moecp]) =0, xeB{S)\S 

and the special principle: 

"Let V be any vector field on ]R" that is tangent at S and also tangent at all surfaces 
displaced by a fixed distance from S (i.e., all surfaces defined as level sets of the 
distance function d to S). Then at points y on the surface divg v{y) — dwv{y)." 

We give a proof of this by appealing to Lemma 4.1 but with the additional assumption of C^- 

regularity of the vector field v if codim S > 2. 

Proof. Given the parametrization 7 of a generic surface patch, the map 

0(01,02) =7(0i)+ No 7(0i)- 

parametrizes ecp^^ (7(^'i)) for fixed 9i and gives us the required coordinate system. For points x 
off the surface the normal Vd to the level-sets of the distance map is given by 



|x-ecp(x)| ' ' ' |N 07(01) -021 

Then the image of the following projection matrix 

Q{x) = l-Vd{x)Vd{x)^ ^ Q(0i,g2) = QocD(0i,02) = J- j^'.^^p 

is tangent to level lines of the Euclidean distance map d. Since v is itself tangent to the latter level 
lines, we have Q ■ c = c or in 0-variables Qv = v with 

i;(0i,02) = 0(01, 02) = Dej7((^i) ■^((^i,^^2) + No7(0i) • J7(0i,02) 

as in Lemma 4.1. Since Vd is also normal to the surface S, the product Qv is 

Qv = De,^ ■ V + QNi] . 

If the codunension is one, the matrix N is an n x 1 matrix and Vd = ±N, hence, the summand 
QNrj cancels out, and this product reduces to Qv = Dg^ 7 ■ v. Consequently, in order for Qv = v 
to hold, T] = must vanish, and Lemma 4.1 yields the result. If the codimension is higher, we 
rewrite the product as 



Qv = De,7-v + N. J - 



020^ N'N 
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The coefficient must satisfy ?/(0i,O) = so that v{x) is tangent to the surface if x G S. In order to 
have additionally Qd = v,r] must also satisfy 

ejAf]=0 where A = := (N o 7(61))^ ■ N o 7(61) . 

We differentiate the new condition 62 At] = with respect to 62 and obtain 

tj'^A + ejADg^rj =0 => // = -A-^Dg^t]^A02 

Since in higher codimension v is assumed to be C^, we can differentiate again which yields 

Dg^>] = -A-^Dg^r]^Ae2-A~^DlT^^Ae2 . 

We apply the trace-operator and get 

trace Dg^rj = - trace Dg^tj'^ - trace (^A~^Dg^;/^A02) ■ 

Finally, since the second derivative Dg^rj is bounded on a small compact neighborhood of (0i,O) 
we have 

trace De,v = 0{\e2\) , 02^0. 
So the second requirement on i] is satisfied and Lemma 4.1 yields the result. □ 

5. Construction of Closest Point Functions from Level-Set Descriptions. Beginning with 
the case of codim S = 1, we present a general construction of closest point functions in the special 
case when the surface is given by a level set (or as an intersection of several level sets.) 

5.1. Codimension One. Let cp : B(S) S denote a smooth scalar level-set function. The 
surface S is the zero-level of cp, which we assume is a proper implicit description of S, i.e., 

V(p{x) 7^ for all x e B{S) . (5.1) 

Thus the normals to level-sets are given by N : B(S) R", N = \7q)/\\7q)\ and the normal 
field of S is N|s : S R". Again B(S) C IR" denotes a band around S and condition (5.1) will 
determine a reasonable band B(S) when starting out with a cp defined on all of M". 

Now closely following [ ], we construct a retraction by solving the following initial value 
problem (IVP): 

^' = -Vcpo^, ^(0,x)=x. (5.2) 

We denote the corresponding family of trajectories by ^(t, x), where t is the ODE-time, while the 
parameter x refers to the initial point. If we start at a point x with cp{x) > and solve forward 
in ODE-time, f is a steepest descent trajectory, while if we start at a point x with cp{x) < and 
solve backward in ODE-time, we will obtain a steepest ascent trajectory heading for the surface. 
The unique intersection of the trajectory f ( ■ , x) with the surface S defines a retraction that maps 
the initial point x to some point on S. The next step is to find the point of intersection which we 
achieve by a suitable transformation of ODE (5.2). We consider the descent case (<p(x) > 0), and 
relate the ODE-time t and the level label A by 

<P(^(t(A),x)) = ^(x)-A. (5.3) 

The implicit function t : [0, cp{x)] — > ]R-|_, A — )• t(A) takes the value t(0) = and has the 
derivative t' = l/\Vcp\'^ o ^. The transformation that we want is ;/(A, x) := ^(t(A),x) because 
(5.3) can be rewritten as ^(;^(A, x)) = cp{x) — A, and thus evaluating rj at X = cp{x) returns the 
corresponding point on the surface S: 



(p{rj{(p{x),x)) =0 ^ ?/((p(x),x) e S. 



CALCULUS ON SURFACES WITH GENERAL CLOSEST POINT FUNCTIONS 



13 




(5.4) 



So far we have considered the descent case but for ascent we end up with the same IVP. Finally, 
we obtain the desired closest point function cp by 



By construction, this function is a differentiable retraction. That it is in fact a closest point function 
is a direct consequence of Theorem 3.7: because of our construction the pre-image cp~^(i/) is 
exactly a trajectory of the ODE (5.4) and this trajectory intersects S orthogonally. 

5.1.1. Remarks. 

• If one is only interested in the construction of a retraction map it suffices to replace in (5.4) 
with a transversal field c, i.e., |c| = 1 and (c, N) > /5 > 0. After a similar transformation the 
retraction is given byr(x) = r]{q){x), x). This approach to retractions is actually the method of 
backward characteristics (see e.g., [ ]): the solution of the PDE 



is V = u o r. So, in a neighborhood of S, v defines an extension of the surface function u and 
is as smooth as u (see the method of characteristics in e.g., [ ]). The idea of extending surface 
functions by solving (5.5) with c = \/cp has been used in other earlier works, see e.g., [ ]. 

• Note that the same ODE as in (5.4) is used in [ ] for the construction of diffeomorphisms. 

• If / : R ^ R is a smooth function with then the closest point function obtained from 
using q> := f o cp in IVP (5.4) is the same as that obtained from the original cp because the 
corresponding ODEs parametrize the same curve through the initial point x. 

• The Euclidean closest point function is a special case our construction: let q) = sd he the 
signed Euclidean distance function. In this case, the IVP for t] reduces to rj' (A, x) = — V sd(x), 
?/(0, x) = X, with a right-hand side independent of A. This has the solution rj{\,x) = x — 
AV sd(x) and the corresponding closest point function is the Euclidean one 



5.2. Higher Codimension. In the case of higher codimension, codim S = m > 2, we as- 
sume that S is the proper intersection of m surfaces of codimension one. Given this, we construct 
the closest point function cp = cp^ o ■ ■ ■ o cpj o cpj as the composition of m closest point func- 
tions, where cp^ retracts onto Si, cpj retracts Si -intrinsically onto Si n S2, cpg retracts Si n S2- 
intrinsically onto Si H S2 H S3, and so on. We demonstrate only the situation where codim S = 2, 
since this is essentially the inductive step for the general case. 

Let again B(S) C M" be a band around S and let : B(S) S,j E {1,2}, be smooth level-set 
functions. Each level-set function q)j shall yield a proper implicit description of a codimension- 
one surface Sj as its zero-level, thus we require: 



The two surfaces Si and S2 are orientable with normal vector fields given by Nj : Sj R", 
Nj = \7 cpj / \\/ cpj\. We assume their intersection S = Si H S2 to be proper, that is, the normal 
vectors have to be linearly independent on S. In fact we assume the linear independence to 
hold on all of B(S) (possibly we have to narrow B(S)). This implies that any two level -sets 
Sj*^ := {x : <j?i(x) = ^1}, Sj^ := {x : cp2{x) = ^2} intersect properly (as long as the intersection is 
non-empty). 

Now, we set up a closest point function cp = cpj o cp^ in two steps. The first step is the same 
as in the previous section: cp^ : B(S) — > B(S) n Si maps onto the subset B(S) n Si of the first 
surface Si by cpj = ?/i((pi(x),x) where rji solves (5.4) with cpi in place of cp. 



cp :B(S)^S, 



cp(x) := y]{cp{x),x) . 



{c{x),Vv) = , v\s , 



(5.5) 



?/(sd(x),x) = X — sd(x)Vsd(x) = ecp(x) . 



Vcpj{x)^0, for all x e B(S) . 
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Fig. 5.1. Example 1 (left): a circle (red) embedded in R^, given as the intersection of a sphere and a plane. Example 2 (right): 
a curve (red) embedded in R^, given as the intersection of a cylinder and a parabola. 

The second step is more interesting: we construct an Si -intrinsic retraction cpj : B(S) n Si 
S. The idea is essentially the same as that of the previous set-up, but now it is Si -intrinsic: we 
consider the IVP 

ri2{Kx) = -——i — r^o;^2(A,x), ri2{0,x) = x , (5.6) 
I 921 

by transforming analogously to the codimension-one case (c.f., (5.4)). Note that Vsj9'2 is given 
by = Pi^92 with = I — N-[NJ. Because we start in Si and Vs^fi is tangential to Si, 

the curve ?;2( ■ , x) is contained in Si. Hence, we obtain an Si -intrinsic closest point function cp2 
by 

cp2 : B(S) n Si — > S , cp2(x) := t]2{cp2{x) , x) . 

The composition cp = cp2 o cpj defines certainly a differentiable retraction cp : B(S) S, and, 
as all the intersections are orthogonal. Theorem 3.7 guarantees cp to be a closest point function. 

5.3. Example 1. We construct a closest point function for a circle S of radius embedded 
in IR'' as the intersection of a sphere and a plane (Figure 5.1 left). The two level set functions are 

(pl{xi,X2,Xo,) =x\ + x\ + xl-l , (p2{xi,X2,X^) = X3 - 1 . 

The equation cpi = yields the unit sphere as the codimension-one surface Si, equation <p2 = 
specifies the second codimension-one surface S2 which is a plane parallel to the XiX2-plane at 
X3 = 1/2, and the circle S = Si H S2 is the intersection. 

First step: we set up a closest point map cpj onto the sphere Si. Let ci = V<pi, the IVP for r]\ 

is 

n[ = -^r^o^^ = -^\^, t^,{Q,x)=xe^^\{Q}. (5.7) 

|V(pi| ^-mr 

The maximal band around the sphere where V(pi 7^ does not vanish is IR'^ \ {0} . The solution 
rji of IVP (5.7) and the corresponding closest point function cpj : ^ \ {0} S\ are 

;?i(A,x) = ^|x|2-A- 1^ ^ cpi(x) = ?7i((pi(x),x) = 1^ . 

(Note that (pi is of the form cpi = / o sd with / strictly monotone. Hence cp^ is the Euclidean 
closest point function, see the remarks 5.1.1.) 
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Fig. 5.2. A circle {red) embedded in R^, given as the intersection of a sphere and a plane. Left: the closest point function 
cp = o cpj maps the black diamonds to the red dots by a two-stage retraction: first cpj maps the black diamonds onto the blue 
crosses on the light blue sphere, and second cpj maps blue crosses onto the red dots on the red circle by following trajectories on the 
sphere. Right: the closest point function cp(x) = cpj o cp.2_[^) first maps the black diamonds onto the blue crosses in the light blue 
plane, then maps the blue crosses onto the red dots on the red circle by following trajectories contained in the plane. 



Second step: we set up an Si-intrinsic closest point map cpj onto the circle S. Let C2 
(P2 = Pi ■ V^2- The projector Pi is given by 



Pi{x)=I- 



X I X^ I 



~1~ — X\X2_ — X\X'^ 

— X1X2 + X3 — X2X3 

^2 1 ^2 



-X1X3 — X2X3 X 



iv<pi(x)r 

For X E Si, i.e., |x| =1, we get then 

Vs^Cplix) = Pi{x) ■ V (p2{x) = (-X1X3, -X2X3, X1+X2) = (^-XiX3, -X2X3, 1-X3^ , 
|Vsi<p2W|2 



1 — X: 



3 • 



The largest subset of Si where Vg^ (p2 does not vanish is the sphere without the poles Si \ {±£3}, 
and so the IVP for rj2 is 



V'2 



° '^^ = ■ ^ ■ ^ ' V2{0,x)=xeSi\ {±es} , (5.8) 



where t]2^j is the j-th component of r]2- The solution rj2 of IVP (5.8) and the corresponding closest 
point function cp2 : Si \ {±£3} S are 



'72 (A, x) 



l-(^3-AF 
1-xl 



X\, 



X2, X3 — A 



cp2(x) = rj2i(p2{x),x) 



^1' 7 X2, 1 



Finally, we compose the closest point function cp : B(S) = IR'' \ {x : xi = X2 = 0} ^ S by 

T 



Cp(x) = Cp2 0Cpi(x) 



1 



Figure 5.2 shows this construction of cp schematically. The maximal band B(S) around S — 
where cp is defined — is M.^ without the X3-axis, since the X3-axis gets retracted by cpi to the 
north-/ south-pole of the sphere Si where cpj is not defined. 
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This first example is intended to highlight the concept of our approach. In fact, in this par- 
ticular case, it is much simpler to first project onto by cp2(^) = (^1/^2/1/2)^, and then retract 
S2-uitruisically (i.e., in the plane) onto S by cpj : S2 \ {0} ^ S, to obtain cp(x) = cpj o cp2(^) 
This approach is also illustrated in Figure 5.2 and, in this particular case, yields the same closest 
point function (and in fact they are both equal to ecp). 

5.4. Example 2. We consider a curve embedded in given as the intersection of a cylinder 
with a parabola (see Figure 5.1 right). The two level set functions are 

fl(xi,X2,X3) = 1 - X2 -X3 , Cp2ixi,X2,X3) = X3- xl . (5.9) 

This time we find the closest point functions numerically using the ODE solver ode45 in 
MATLAB. The first closest point function cp — which first maps onto the cylinder — is obtained by 
solving the two ODEs below 

^1 = - |^^^|2 ° Vi ' f/i(0,x)=x, ^ z := rii{cpi{x),x) , 



|V<p 



ri2 = -j^^^^°V2, ?72(0,z)=z, cp{x) ■.= r]2{q>2{z),z) , 

in the given order. In the same way we obtain the second function cp — which maps first onto the 
parabola — ^by interchanging the roles of cpi and cp2 and numerically solving the ODEs 

V2 = - ° V2 , //2(0,x)=x, =^ z--r]2{(f\(x),x), 

\y (f2\ 

Vi = - i^^VyiV ° ' ^ cp{x) ■.= rii{cpi{z),z). 

The retraction stages of the resulting closest point functions cp and cp are visualized in Figure 5.3. 

In contrast with our first example, here the two closest point functions cp and cp are different. 
We define a 50 x 50 x 50 Cartesian grid G on a reference box R = [—1.25, 1.25] x [—1.25, 1.25] x 
[—0.25, 1.25] containing S. We use a subset of these points as a narrow band of grid points sur- 
rounding Sr Comparing the values of the distance function q) = {q>^ + (P2) ^ at the closest points 

max (p( cp(x) ) = 4.4893 ■ 10~^^ 

xeGnB{s) 

max cp( cp(x) ) = 4.4758 ■ 10"^^ 

xeGnB{s) 

we can see that both numerical functions cp and cp are very accurate. We also note they are in- 
deed different mappings because max^ggpg^g) ] cp(x) — cp(x) ] = 1.7095 ■ 10^'^'' (see also Table 6.2 
where they are clearly distinct from ecp). 

6. The Closest Point Method with Non-Euclidean Closest Point Functions. In this section 
we demonstrate that the explicit Closest Point Method based on Euler time-stepping still works 
when replacing the Euclidean closest point function with non-Euclidean ones. 

Given an evolution equation on a smooth closed surface S as 

dtu — As{t,y,u) — , u{0,y) = UQ{y) , yES 

where A5 is a linear or nonlinear surface-spatial differential operator on S, following °], the 
semi-discrete explicit Closest Point Method based on Euler time-stepping with time step t is 

Initialization: vq = uqo cp 

Evolve step: = + 'rA(f„,x, , x E B{S) 

Extension step: ^n+i = ^n+i ° cp 



^For this particular example we define a band by B(S) — {(xi,X2, -T3) : cp < 0.125} where we use (p = + "Pi)^ ^ 
lieu of Euclidean distance. Our banded grid is G fl B(S) and contains 1660 grid points. It contains more points than are 
strictly necessary: see [ , Appendix A] for an approach to banded grids for the Closest Point Method. 
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Fig. 5.3. A curve {red) embedded in V? , given as illustrated in Figure 5.1 (right). Left: the closest point function cp = 
o cpj maps the black diamonds to the red dots by a two-stage retraction: first cpj maps the black diamonds onto the blue crosses 
on the light blue cylinder, and second cpj maps the blue crosses onto the red dots on the red curve by following trajectories on the 
cylinder surface. Right: the closest point function cp = cpj o cp,^ first maps the black diamonds onto the blue crosses on the light 
blue parabola by cpj, then cpj maps the blue crosses onto the red dots on the red curve by following trajectories on the parabola. 



where A is a spatial operator on B(S) which is defined from Ag via the closest point calculus and 
hence 

A{t,x,uoc:p)\x=y = As{t,y,u) , yeS. 

For the fully discrete Closest Point Method which is also discrete with respect to the spatial vari- 
able X E B(S) we replace A with a discretization of it, while the extension tf„+i o cp(x) is replaced 
with interpolation of the discrete w„^i in a neighborhood of cp(x) because the point cp(x) is not 
a grid point in general [ ] . 

6.1. Advection Equation. We consider an advection problem with a constant unit speed on 
the curve S shown in Figure 5.1 (right): 

3fM + divg (m ■ T(i/)) = , m(0,i/)=i/3, y^S, 

V<pi(y) X v<p2(y) (6.1) 



with T{y) 



\Vcpi{y) X V(p2{y)\ 



The surface-spatial operator is As{y, u) = — divg [u ■ T{y)), while the operator used in the Clos- 
est Point Method is 

A{x,v) = — div{v ■ T o cp{x)) , v = uocp, x e B(S) 

according to the Divergence Principle 3.5. 

We solved this advection problem with three different closest point functions — ^namely cp, 
cp from Example 2, and the Euclidean closest point function ecp — and with three different mesh 
sizes on the reference box _R = [—1.25,1.25] x [—1.25,1.25] x [—0.25,1.25]. The ecp was com- 
puted with a numerical optimization procedure using Newton's method. The evolve step uses 
the first order accurate Lax-Friedrichs scheme (time step t = 0.95/z in accordance with the 
CFL-condition) and the extension step is performed with WENO interpolation [ ] (based on 
tri-quadratic interpolation). The solution is advected until time t = 1. 

We compare to a highly accurate solution obtained from parametrizing the problem. By using 
the parametrization 7 : [0,2n) — >■ S,'y{6) = (cos(0),srn(0)-\/l + cos(0)^,cos((>)^)^ the advection 
problem (6.1) is equivalent to 

l-y' {0)\dtu - dgu = , u{0,e) =cos{9f , u{t,2n) = u{t,0) , (6.2) 
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Table 6.1 

Errors measured in loo-norm at stop time t = 1/or Closest Point Method approximation of the advection problem (6.1) on the 
curve shown in Figure 5.1 (right) using various closest point functions. There is no significant difference regarding the choice of the 
closest point function. 



h 


cp 


cp 


ecp 


0.0125 
0.00625 
0.003125 


9.1630e-03 
4.6182e-03 
2.3262e-03 


9.1592e-03 
4.6172e-03 
2.3260e-03 


9.0615e-03 
4.5544e-03 
2.2908e-03 



where u{t,6) = u{t,y{9)). The formal solution of the latter is 

e 

u{t, 6) = cos (s-^ {t + s{9))'j , where s{9) = J \Y {ix)\ da. (6.3) 



is the arc-length function. We used Chebfun [ ] within MATLAB to find a highly accurate ap- 
proximation to s{6) and thereafter a Newton iteration to approximate the value s^^ {t + s{6)). 

Table 6.1 shows the error in the results measured in /t»-norm at f = 1. We see there is no 
significant difference regarding the choice of the closest point fimction and that the error is of 
order 0{h) as expected from the Lax-Friedrichs scheme. 

6.2. Diffusion Equation. Here, we consider the diffusion equation on the curve S shown in 
Figure 5.1 (right): 

dtU-Asu=0, uiO,y) = ^^P(fr) ^ ygg^ (6.4) 

50 

The surface-spatial operator is now Ag (y, m) = A5M. 

Here again, we compare to a highly accurate solution obtained from parametrizing the prob- 
lem. By u{t,6) = u{t,j{6)), where 7 is the same parametrization that we used in (6.2), the 
diffusion equation (6.4) transforms to 

The formal solution of the latter (with frequency parameter co = 2n/s{2Tz) and arc-length func- 
tion s(0)) is given by 

m=-oo S[ZTl) J DU 

We again used Chebfun [ ] to obtain highly accurate approximations to c„, and s{6). Moreover, 
we restrict the summation to — M < m < M where M is chosen such that the bound e~'^^*^^'co is 
lower than machine accuracy. 

In Section 4, we have discussed three different ways to deal with the Laplace-Beltrami oper- 
ator by the closest point calculus which in turn give rise to three different operators A. Now, we 
solve the diffusion equation (6.4) using the Closest Point Method with each of these (in the evolve 
step) and with each of our three different closest point functions ecp, and cp, cp from Example 2. 
The extension step is performed with tri-cubic interpolation. The time step is chosen as t = 0.2h^ 
for numerical stability. The errors are measured at time t = 0.1 in /00-norm. 
1. Writing the Laplace-Beltrami operator as A5M = div5(VsM), the simplest possibility is 

A(z;) = div(Vz^) = Az; . (6.5) 

But, recall that this need not work if the closest point function does not satisfy the requirement 
of Theorem 4.3. A is discretized by the usual -accurate 7-point stencil. Table 6.2 shows 
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Table 6.2 

Errors measured in Irx-norm at stop time t = 0.1 for Closest Point Method approximation of the diffusion equation (6.4) on the 
curve shown in Figure 5.1 (right) using various closest point functions. The spatial operator on the embedding space is A{v) = Av 
as in (6.5). We observe convergence when using ecp as expected. But, here, cp and cp do not yield a convergent method as they (with 
hindsight) do not satisfy the requirement of Theorem 4.3. 



h 


cp 


cp 


ecp 


0.0125 
0.00625 
0.003125 


0.017394 
0.017499 
0.017526 


0.017394 
0.017499 
0.017526 


4.717879e-04 
1.173944e-04 
2.927441e-05 



Table 6.3 

Errors measured in l^o-norm at stop time t = 0.1 for Closest Point Method approximation of the diffusion equation (6.4) 
on the curve shown in Figure 5.1 (right) using various closest point functions. The spatial operator on the embedding space is 
A(x, v) = div(P o cp(a:) • \/v) as in (6.6). There is no significant difference regarding the choice of the closest point function. 



h 


cp 


cp 


ecp 


0.0125 
0.00625 
0.003125 


4.805473e-05 
1.198697e-05 
2.975015e-06 


4.856909e-05 
1.211579e-05 
3.008249e-06 


4.793135e-05 
1.195396e-05 
2.968171e-06 



the errors. We observe convergence at the expected rate 0{h^) when using the Euclidean 
closest point function ecp, while cp and cp do not yield a convergent method. With ecp, as 
opposed to cp and cp, we can be sure that = Ac, since ecp satisfies the requirement of 
Theorem 4.3. 

2. Now, we rewrite the Laplace-Beltrami operator as Agw = divg (PVsm), where the projector P 
is given by the outer product P{y) = T{y) ■ T{yY and T is the tangent field as in (6.1). Based 
on this formulation we obtain the following operator 




and Theorem 4.2 ensures that the operators coincide on the surface, i.e., A\g = As for all 
closest point functions in accordance with Definition 3.1. A is discretized by replacing the par- 
tial differential operators with corresponding C)(/i'^) -accurate central difference operators. 
Table 6.3 shows the errors. We observe convergence at the expected rate 0{h?-) for all three 
closest point functions and there is no significant difference regarding the choice of the closest 
point function. Note that this approach does require that we know the tangent field. 
3. Finally, we work with re-extensions and appeal to the Principles 3.4 and 3.5: 

3 

A{v) = div(Vz; o cp) = ^ d^^ [{d^^v) o cp] . (6.7) 
1=1 

The discretization of A involves two steps: the partial differential operators 9x, are replaced 
with corresponding central difference operators, the re-extensions o cp are replaced with 

tri-cubic interpolation. Table 6.4 shows the errors. We observe convergence at the expected 
rate 0{h^) for all three closest point functions. There is no significant difference regarding the 
choice of the closest point function. 

7. Conclusions. We presented a closest point calculus which makes use of closest point func- 
tions. This calculus forms the basis of a general Closest Point Method along the lines of [ 18]. The 
original Closest Point Method of [ ] inspired the present work: after becoming aware of the 
key property D ecp(i/) = P{y) which accounts for the gradient principle of the original method, 
we turned it into a definition of a general class of closest point functions namely Definition 3.1. 
Here, we characterize such closest point functions (see Theorems 3.3 and 3.7) and prove that this 
class of closest point functions yields the desired closest point calculus (see Theorem 3.2 and the 
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Table 6.4 

Errors measured in loo-norm at stop time t = 0.1 for Closest Point Method approximation of the diffusion equation (6.4) 
on the curve shown in Figure 5.1 (right) using various closest point functions. The spatial operator on the embedding space is 
A{v) = div(Vi; o cp) as in (6.7). There is no significant difference regarding the choice of the closest point function. 



h 


cp 


cp 


ecp 


0.0125 
0.00625 
0.003125 


1.508222e-04 
3.754836e-05 
9.347565e-06 


1. 50341 le-04 
3.743075e-05 
9.319686e-06 


1.492513e-04 
3.715598e-05 
9.251365e-06 



Gradient 3.4 and Divergence Principles 3.5 which follow). The closest point calculus is also suf- 
ficient to tackle higher order surface differential operators by appropriate combinations of the 
principles, however for all surface-intrinsic diffusion operators the calculus can be simplified to 
use less closest point extensions (see Theorems 4.2 and 4.3). Finally, the basic principles of the 
original Closest Point Method are now contained and proven in our general framework. 

In addition to the framework, we describe a construction of closest point functions given a 
level-set description of the surface and show examples. This demonstrates that there are inter- 
esting closest point functions besides the Euclidean one. Furthermore, we demonstrate on two 
examples (the advection equation and the diffusion equation on a closed space curve) that the 
Closest Point Method combined with non-Euclidean closest point functions works as expected. 
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